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Abstract 

Multi-peaked localized stationary solutions of the discrete nonlinear Schrodinger 
(DNLS) equation are presented in one (ID) and two (2D) dimensions. These are 
excited states of the discrete spectrum and correspond to multi-breather solutions. A 
simple, very fast, and efficient numerical method, suggested by Aubry, has been used 
for their calculation. The method involves no diagonalization, but just iterations of 
a map, starting from trivial solutions of the anti-continuous limit. Approximate 
analytical expressions are presented and compared with the numerical results. The 
linear stability of the calculated stationary states is discussed and the structure of 
the linear stability spectrum is analytically obtained for relatively large values of 
nonlinearity. 
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1 Introduction 



The DNLS equation, Eq. (1), has been extensively used as a generic model 
for studying nonlinear effects (breathers, for example) in a discrete system 
[1,2,3,4,5,6,7,8,9,10,11,12]. In addition, specific applications have been pro- 
posed for the description of: i) local intramolecular stretching vibrations in 
symmetric polyatomic molecules [13], ii) arrays of coupled nonlinear optical 
waveguides [14], in) interacting electron-lattice models [15,16] (or, equiva- 
lently, intramolecular excitation-phonon coupled systems [17,18,19]) in solid 
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state physics, where localized solutions of DNLS correspond to polarons (vi- 
brational polarons, respectively), and recently iv) Bose-Einstein condensates 
[20,21]. 

In DNLS the evolution of a complex probability amplitude \I/ n at the site n of 
a d- dimensional lattice is given by 

iT^^E^ + xWn. (!) 



where V represents nearest neighbor coupling, \ is the strength of nonlinearity, 
and the sum over 5 contains the nearest neighbors of the lattice site n. For 
example, in ID is J2s ^n+s = ^n+i + ^n-i, while in a 2D square lattice, if 
each lattice site n is represented by a pair (n x ,n y ), is J2s ^n+s = ^n x +i,n y + 

^n x -l,n y + ^n x ,n y +l + ^n x ,n y -l- 

The DNLS equation is derived from the Hamiltonian (assuming an infinite 
lattice, or periodic boundary conditions) 

^-^EE*:^ + fEI^I 4 - ( 2 ) 

n S n 



There exist two conserved quantities during the DNLS dynamics: the Hamil- 
tonian (2) and the norm 

N = T,\^n\ 2 - (3) 



A trivial transformation (\I/ n — > ty n / \/~N) connects the solutions of DNLS for 
an arbitrary norm N with those normalized to unity, through a rescalling of 
the nonlinearity (x ~ * X ' N). Therefore, any solution $ n of Eq.(l) with norm 
N is obtained from a solution \P„ with norm 1, through 

^ n (t;V,x) = VN-^n(t;V,xN) (4) 



Furthermore, the magnitude of the hopping integral V can be considered unity 
by appropriate rescalling of time and x- Changing sign in V is equivalent to 
the transformation * n -> (-l) n *„ (or V nx ,n y -> (-^-) nx+ny ^n x ,n y in 2D, etc.). 
Therefore, in the following we consider N — 1 and V = 1 without loss of the 
generality, while the nonlinearity x remains the only free parameter of the 
system. 
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1.1 Stationary solutions of DNLS 



The stationary states of DNLS are characterized by a simple harmonic evolu- 
tion with frequency uj: 



^n=i>n- e-™. (5) 



The time independent amplitudes ip n satisfy the nonlinear eigenvalue problem 

^i>n = -J2^n+S + X\^n\ 2 ^n (6) 
<5 



(since we consider V — 1). The energy E s of a stationary state is related to 
its frequency uj through 

£ s =u;-f£ivg 4 . (7) 



There are two kinds of stationary solutions of Eq. (6): extended (like Bloch 
states and standing waves [22]) and localized states. The Bloch states, ip^ ~ 
exp(ign), form a band of frequencies or energies from —2dV to 2dV (for an in- 
finite c?-dimensional lattice). Localized states have discrete frequencies outside 
of the Bloch band. 

The most obvious and well studied localized state is the single-peaked solu- 
tion, which has its maximum amplitude on one lattice site. This state has 
extreme (minimal for negative x an d maximal for positive x) energy and 
frequency, compared to the other localized states. However, in general there 
are infinite (in an infinite lattice) multi-peaked stationary solutions of DNLS 
(which, for example, can be continued from trivial multi-peaked solutions of 
the ant i- continuous limit), resulting in a very rich discrete spectrum with 
many quasi-degenerate levels (see below). In spite of this complexity, a part of 
the spectrum can be understood by classifying the stationary states from the 
anti-continuous limit. For the general concept of the anti-continuous, or anti- 
integrable, limit in nonlinear lattices see Ref. [23]. Here, multi-peaked and real 
stationary solutions of high symmetry, which have direct counterparts at the 
anti-continuous limit, are presented. Their location on the DNLS spectrum, 
their stability, and approximate analytical expressions are discussed. 
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2 Aubry's method for the numerical calculation of localized sta- 
tionary states of DNLS 



Localized stationary states can be obtained as attractors of the map: 

4, — ► W= sgnx • (8) 



In this equation, tp = (ipi, . . . , ipi), where L is the total number of lattice sites, 
sgnx denotes the sign of x, ^{4>} is defined through the right- hand- side of the 
stationary equation (6), with its n-th component given by 

H^}n = ~Y1 ^™+<5 + X WVn, (9) 
<5 



and ||N{V'}|| represents its norm yJ2n=i l^{^}n| 2 - 

Starting the iterations from appropriate initial states, obtained through trivial 
solutions of the anti-continuous limit, i.e. 

^ r=0) = <Wo, or Vi r=0) = -40Sn.no ±tf B>Bl ), etc., (10) 



depending on the desired localized stationary solution (single-peaked, double- 
peaked with interpeak separation \no — ni\, etc.), and iteratively applying the 
map (8) 

^(r+1) = SQny . ^ K ( U ) 



the procedure can rapidly converge to the corresponding localized state. 

Up to now this method has been successfully applied for calculating single- 
peaked stationary states of DNLS in ID, 2D, and 3D [16], as well as for 
the single-peaked ground states in other similar systems [24,25,26,27,28]. The 
method has been invented by S. Aubry for finding polarons in the adiabatic 
Holstein model [29,30,31], a problem which reduces to the stationary solutions 
of DNLS. 

In the following, after briefly recalling some results obtained in Ref [16] re- 
garding the stable single-peaked stationary states of DNLS (section 3), this 
method is applied for calculating multi-peaked stationary solutions in ID (sec- 
tion 4) and 2D (section 5). Stationary states are presented in sections 4 and 
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5 for negative values of x\ the obtained solutions in this case can be directly 
transformed to the corresponding ones at — x (i-e. at positive nonlinearities) 
by changing the sign of the frequency (ou — > —uo) and energy (E s — > —E s ), 
and making the transformation ^ n — >• (-l) n ip n (or ^ nx . riy -> (-l)"^™^^ 
in 2D, etc.) in the wavefunction. However, analytical results and the general 
discussion concern both signs of X- 



3 Single-peaked (SP) stationary states 

For negative (positive) nonlinearity the single-peaked solution of DNLS corre- 
sponds to the lowest (highest) frequency stationary state. In ID there is always 
a SP state with extreme frequency and energy, for any nonzero value of x- As 
X is approaching zero from negative (positive) values, the frequency and the 
energy of the SP solution tends to the bottom (top) of the Bloch band. The 
branch of SP solutions, as obtained by varying x, has two well-known limits: 
for large \x\ (in the anti-continuous limit, where the first term of the right- 
hand-side of Eq. (6) can be neglected) tends to a single-site localized state 
(ip n = 5n,n ), while for Ixl — > tends to a solution obtained by the soliton of 
the continuous nonlinear Schrodinger equation (see Eq. (14) below). 

The picture is qualitatively different in 2D and 3D. In these cases there is a 
critical value of nonlinearity xi > 0, such that for \x\ < Xi does not exist a 
single-peaked, or any other localized, stationary state. At \x\ = Xi a P a h °f SP 
states appears, through a saddle-node bifurcation; a narrow stable solution of 
high amplitude and an unstable one of relatively large extent and small am- 
plitude. Due to the simultaneous existence of two SP states and the proximity 
of the unstable one with the extended Bloch states of the band edges, it is not 
paid any attention to the unstable solution and from now on we exclusively re- 
fer to the stable one wherever a single-peaked state is mentioned in 2D or 3D. 
A second nonlinearity threshold xi > Xi exists, such that for xi < \x\ < Xi 
the SP state has extreme frequency, but not extreme energy, since its energy 
E s lies inside the Bloch band. Only for \x\ > X2 the single-peaked stationary 
solution provides an extreme of the energy (i.e. it is the ground state for neg- 
ative x)- hi 2D the values of these thresholds are xi ~ 5.701 and X2 ~ 6.679, 
while in 3D are xi ~ 7.852 and X2 ~ 10.816 [16] 1 . 

Analytical approximate expressions have been presented in Ref. [16], which ac- 
curately describe the SP stationary states. In particular, for the whole branch 
of SP solutions in 2D and 3D, as well as for values of |x| larger than about 3 



1 Note that the nonlinear parameter x used in this work corresponds to —k 2 of 
Ref. [16]. 
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in ID, the exponentially decaying function 



t'Zn y ,«, = (i^r ■ <'-' + '-'+'-', With C = -I - ^ (12) 



(where one has to disregard the index n z in 2D and both n z and n y in ID), 
can be used to describe the exact SP solution. The corresponding frequencies 
and energies are given by 



u = x-^— 3 and E a = $ + — + 1 3 J . 13 



The expression (12) is derived from a variational method, by employing a 
perturbative expansion of ( in powers of 1/x in the condition providing the 
minima of the variational energy The next non-zero correction of ( in Eq.(12) 
is of the order of 1/x 5 - 

The above analytical results describe accurately (the larger the \x\ the better 
the approximation) the SP solutions of DNLS in all cases, except for relatively 
small values of |%| in ID. Then a smooth transition occurs for \x\ in the region 
2.5 — 3, from the above expressions to the static soliton of the continuous 
nonlinear Schrodinger equation. Therefore, for \x\ smaller than about 2.5 the 
SP solutions in ID can be approximated by 




with corresponding frequency and energy 



u = sgn X ■ (2 + ^ J and E s = sgn X • ( 2 + g ) . (15) 



Note that the {—sgnx) n term in (14) provides the alteration of signs in suc- 
cessive lattice sites, which characterizes the solution at positive x- I n Eq. (12) 
this is obtained through the negative sign of (. 
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Fig. 1. Frequencies of single- double- and triple-peaked stationary solutions of DNLS 
in ID (points). Dashed lines show analytical expressions obtained for large values 
of \x\- The horizontal line at u) = —2 indicates the lower edge of the band of Bloch 
stationary states, which extends from —2 to 2. The spectrum is antisymmetric on 

x; ^(-x) = -<*>(x)- 

4 Multi-peaked solutions in ID 



4-1 Frequency spectrum 



As it has been already mentioned in the introduction, the frequency (or energy) 
spectrum of DNLS in ID comprises many discrete levels, corresponding to 
single-peaked and multi-peaked localized stationary solutions. 

These discrete levels have well determined limits in the anti-continuous regime. 
In this limit the discrete spectrum is given by 

u = ^-, E s = ^-, where M= 1,2,3,... . (16) 



Each (highly degenerate) level of Eq.(16) corresponds to any M-peaked sta- 
tionary state, where all the M peaks (at arbitrary sites) have the same norm 
l/y/M (and arbitrary complex phases). Such stationary solutions can be con- 
tinued away from the anti-continuous limit [30] up to some value of de- 
pending on the particular state. 

A stationary solution at a given value of x is classified as a single- double- or, in 
general, M-peaked state, depending on how many sites (one, two, or in general 
M, respectively) are occupied at the anti-continuous limit of the branch in 
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which this state belongs. Such a classification is always possible for stationary 
states of high symmetry. Moreover, it can be used for each stationary solution, 
providing a complete description of the discrete spectrum, up to u ~ —5.45 
[32]. 

Fig. 1 shows a part of the frequency spectrum which contains contributions 
from many double-peaked solutions, as well as few branches (just shown for 
indication) of triple-peaked states. Analytical expressions, obtained at large 
values of ]x| and shown by dashed lines, describe well these branches (the next 
corrections are of the order of inverse powers of x, see Appendix). Bifurcations 
lead to the disappearance of some branches (or merging with other branches) 
by decreasing the strength of nonlinearity. 

What appears as a middle branch of the double-peaked (DP) solutions in Fig. 1 
actually consists of many branches of closely spaced levels (corresponding to 
DP solutions with any interpeak separation larger than one lattice constant), 
which merge to the level of Eq.(16) for M = 2 as \x\ increases. The lower 
and upper branch of the DP solutions are single branches corresponding to 
symmetric and antisymmetric, respectively, double-peaked states with their 
peaks at neighboring lattice sites (examples are shown at the upper left plots 
of Figs. 2 and 3, respectively). The interpeak separation of a DP state is 
determined by the distance S of the sites where the peaks appear. Symmetric 
and antisymmetric DP solutions in ID, with various interpeak distances, are 
presented in the following subsection. 



4-2 Double-peaked symmetric and antisymmetric solutions 



4-2.1 Stationary states 

In order to find the branches of symmetric and antisymmetric DP solutions, 
the method described in section 2 is applied using as initial states 

^ r=0) = ^ 6n ' ni + 6n ' m) ^™ =0) = ^ 6n ' ni ~~ 6n ' n2 ^ (17) 

respectively. The distance \n 2 — ni\ determines the interpeak separation S of 
the corresponding solution. 

Stationary solutions with interpeak separations 1, 2, 10, and 20 sites are shown 
for different values of x m Fig- 2 for symmetric and in Fig. 3 for antisymmetric 
states. As it is expected, by decreasing \x\ the solutions spread more and more, 
until a value of x where a bifurcation occurs resulting in the disappearance 
of the corresponding branch. The smaller values of \x\ shown at each plot 
of Fig. 2 (apart from the case of S — 1, displayed in Fig. 2a) and Fig. 3 
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Fig. 2. Double-peaked symmetric stationary solutions of DNLS in ID (points) for 
different values of the nonlinearity \ an d various interpeak separations 5: (a) 5 = 1 
lattice site, (b) 5 = 2 sites, (c) 5 = 10 sites, and (d) 5 = 20 sites. Dashed lines 
show analytical approximations of the solutions using Eq. (18) for the discrete cases 
where |x| > 6, Eq. (21) for x = —5.2 in (c) and x = ~ 4 in (d), and Eq. (14) centered 
in the middle between the sites of maximum amplitude for \ = —0.5 up to —5 in 
(a) (see text). 

are close to the bifurcation point and therefore represent one of the most 
extended solutions of the corresponding branch. The case of symmetric states 
with interpeak distance equal to one lattice constant (Fig. 2a) is exceptional in 
this respect, since the corresponding branch survives until the limit \x\ — > 0, as 
it merges to the SP branch for sufficiently small values of As the continuous 
limit is approached, the DP symmetric state with S = 1 becomes practically 
indistinguishable from the SP state, i.e. the solution given by Eq. (14). 

The symmetric (antisymmetric) states with S = 1 provide the lower (upper) 
single branch of the DP solutions of Fig. 1, with frequency equal to uj — | — 1 
(u = f + 1), for relatively large |%|. All the other symmetric and antisymmetric 
DP solutions with interpeak separations S > 1 are congested in the middle 
branch (with frequency equal to uj = | for large In particular, in Fig. 1 
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Fig. 3. Double-peaked antisymmetric stationary solutions of DNLS in ID (points) 
for different values of the nonlinearity x an d various interpeak separations S: (a) 
5 = 1 lattice site, (b) S = 2 sites, (c) S = 10 sites, and (d) S = 20 sites. Dashed 
lines show analytical approximations of the solutions using Eq. (18) in the discrete 
cases where |x| > 6 and Eq. (21) in the two more extended cases of (c) and (d). 

the middle branch contains solutions for S = 2, 3, 4, 5, 10, 20, and 200. 
As the interpeak separation S increases, the two peaks start to not overlap 
much, even for small values of |%|, and in this case the corresponding branch 
survives longer (i.e. persists closer to % = 0). These branches, as \x\ decreases, 
eventually deviate from the line u — and for even smaller values of |x| 
they can be described by double-peaked solutions of the continuous nonlinear 
Schrodinger equation. 

Approximate analytical expressions, which describe the DP states in the most 
of the cases, can be obtained by appropriate superpositions of the single- 
peaked states (12) and (14). One has to take into account that when the two 
peaks do not significantly overlap the norm of each peak is about half of the 
total norm. Then from Eq. (4) follows that the two individual wavefunctions 
superimposed in a double-peaked solution should be provided by Eqs. (12) 
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or (14) corresponding to |. As a result a symmetric or antisymmetric DP 
solution with interpeak separation S, where the first peak is located at the 
site ri\ and the second at 112 = n\ + S (S is assumed to be positive), can be 
approximated by 



tf P = . 1 J \~S~ 2 (C |n - ni1 ± {sgn X ) S &- n ^) , (18) 
^2(1 ± P) V 1 + C 2 v y 



(i + 5)C 5 -(5-i)C 5+2 



or 



where P = {-sgnx) TTC ^ 

md C = -^2-W = "-?' (2 °» 



,pp _ {sgnxT ni \x\ I 1 ± (~g^x) 5 \ , . 

^2(1 ±P) V 16 ^cosh^ii cosh J ' 1 J 

where P = (-sgny) 3 § -^r. (22) 

V ; sinh^ V ; 



The former (latter) solution can be used for a DP state with discrete (rather ex- 
tended) peaks, valid for relatively large (small) values of \x\ - Roughly speaking 
the transition from one form to the other occurs for \x\ in the region 5 — 6. The 
normalization factors P in Eqs. (19) and (22) result from the non-orthogonality 
of the superimposed wavef unctions, i.e. P = J2 n Vn^'^d )V ; f P ' ni+S ''(f )> where 
^SPM^x) denotes the SP solution centered at the site m and calculated for 
the value | of the nonlinearity parameter. The plus (minus) signs in these 
approximate solutions correspond to symmetric (antisymmetric) states, ex- 
cept for the case of positive x an d odd S, where they give the antisymmetric 
(symmetric) DP state. 

Eq. (21) is not applicable for small values of \x\ in the exceptional case of 
symmetric states with S = 1, due to the significant overlap of the two non- 
distinguished peaks (Fig. 2a for \x\ approximately less than 5-6). In this case 
the corresponding DP solution, whose branch is merging to the single-peaked 
branch, can be well described by Eq. (14) centered in the middle between the 
consecutive sites of the two peaks. The analytical approximations discussed 
above have been plotted in Figs. 2 and 3 (dashed lines) along with the nu- 
merical solutions (points) for comparison. In particular, Eq. (14) has been 
used, as just explained, in Fig. 2a for x — —0.5, — 1, —2, and —5, Eq. (21) 
has been used only for the smallest values of |x| in Figs. 2c, 2d, 3c, and 3d 
(where \x\ < 6), and Eq. (18) has been plotted in all the other cases. The 
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more distinguished the two peaks, i.e. the higher the S, or the higher the \x\ 
even for smaller values of S, the better the approximations (18) and (21) are. 

4-2.2 Linear stability 

Symmetric and antisymmetric DP states show qualitatively different behavior 
regarding their stability. A detailed investigation of the stability eigenvalues 
is presented below for negative values of x- From this, the behavior at pos- 
itive x can be obtained as follows: for even S the situation, regarding the 
stability eigenvalues, is exactly the same as that of x < 0, while for odd S 
the symmetric (antisymmetric) states behave exactly like the antisymmetric 
(symmetric) states of x < 0. Note here that the transformation (— l) n ip n , con- 
necting stationary solutions at opposite values of Xi turns a symmetric DP 
state to antisymmetric and vis versa when S is odd. Therefore, what is men- 
tioned below for x < 0, also holds as it is for x > when S is even, while it 
valids after interchanging roles between symmetric and antisymmetric states 
when S is odd. 

The eigenvalues of the linear stability problem are always obtained as pairs 
of opposite sign and there always exists a pair of eigenvalues at zero. In our 
approach (see Appendix), if there is an eigenvalue with non-zero imaginary 
part, then the stationary solution is unstable. All relevant discussions in this 
and the following section will refer to those eigenvalues with non-negative real 
part, without usually mentioning the pinned pair of eigenvalues at zero. The 
spectrum of eigenvalues is always symmetric with respect to the imaginary 
axis. The complete structure of the linear stability spectrum close to the anti- 
continuous limit is analytically derived in the Appendix. 

For negative x, the nonlinear term of DNLS has the same sign as the tun- 
neling term. Then it is known that the symmetric DP states are unstable 
[33,34], in contrary to the case of the antisymmetric ones that may be linearly 
stable [4,35]. In fact, the latter are linearly stable in the larger part of the 
corresponding branch. When this branch becomes unstable, by decreasing |x|, 
soon it disappears. The larger the interpeak distance, the smaller the |x|'s 
at which the branch turns unstable and disappears. The scenario for devel- 
opment of instability and disappearance seems to roughly be as follows. For 
large |x|, linear stability analysis provides, a discrete eigenvalue outside of the 
band. For relatively large \x\ the band extends from ^ — 2 to ^ + 2, apart 
from the case of symmetric (antisymmetric) states with interpeak separation 
S = 1, where it extends from ^ — 1 to ^ + 3 (from ^ — 3 to ^ + 1). As \x\ 
decreases, the band moves towards zero, while the discrete eigenvalue keeps 
off zero for S > 2, remains almost constant for S — 2, and goes to zero (but 
slower than the band) for 5=1. After their collision (or the collision of the 
discrete eigenvalue with another eigenvalue splitted off the band) the insta- 
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Table 1 

Nonlinearity strength \ f° r the development of instability and disappearance of 
antisymmetric and the disappearance of symmetric double-peaked stationary solu- 
tions of DNLS at different interpeak separations. Values in parentheses in the second 
column show rough analytical estimates using Eq. (24). 



Interpeak 
separation 


Regime of x at which 
the antisymmetric branch 
becomes unstable 


Regime of x at which 
the antisymmetric 
branch disappears 


Regime of x at which 
the symmetric 
branch disappears 


S= 1 


[-19.73,-19.72] (-18) 


[-9.58,-9.57] 




S = 2 


[-8.8,-8.7] (-8) 


[-6.54,-6.53] 


[-9.63,-9.62] 


S = 3 


[-7.1,-7.0] (-6.3) 


[-5.9,-5.8] 


[-7.9,-7.8] 


S = A 


[-6.5,-6.4] (-5.5) 


[-5.7,-5.6] 


[-7.0,-6.9] 


S = 5 


[-6.2,-6.1] (-5.0) 


[-5.5,-5.4] 


[-6.5,-6.4] 


S= 10 


[-5.0,-4.9] (-4.2) 


[-4.9, -4.8] 


[-5.2,-5.1] 


S = 20 


[-3.90,-3.89] (-4.0) 


[-3.90, -3.89] 


[-4.0,-3.9] 



bility develops. As \x\ decreases more, the band is approaching zero and the 
branch disappears when an eigenvalue splitted from the band collides with 
the pinned eigenvalue at zero. Table 1 shows for antisymmetric DP solutions 
of various interpeak separations (first column) the nonlinearity regime where 
the corresponding branch becomes unstable (second column; for larger values 
of \x\ the solutions are linearly stable) and the nonlinearity regime where the 
branch disappears (third column). 

Regarding the unstable symmetric DP states, linear stability analysis reveals 
that there is always one pair of purely imaginary unstable eigenvalues. If the 
magnitude of instability is denoted by X u (i.e. the unstable eigenvalues are 
±i\ u ), then for fixed value of Xi ^« decreases as the interpeak separation 
increases. For fixed interpeak distance S, but larger than two lattice sites, \ u 
decreases by increasing |x|. For S = 2, \ u tends to a constant value for large 
|x|, while for S — 1, X u decreases with decreasing \x\, in accordance with the 
fact that this branch merges to the stable single-peaked branch for \x\ — > 0. 
This behavior is shown in Fig. 4. In some cases, especially for large interpeak 
separations, the magnitude of instability is so small that for any practical 
purpose the corresponding solution can be considered as quasi-stable [33]. 
The log-log plot of \ u as a function of \x\ in Fig. 4 demonstrates a power-law 
dependence X u ~ \x\ a at large values of \x\- Relating the stability eigenvalues 
with the energy spectrum of a tight-binding problem in the presence of a deep 
and narrow double- well potential (see Appendix), one obtains 

X U = 2 S / 2 Ixl 1 "!, for | X |»1. (23) 
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Nonlinearity |%| 

Fig. 4. Magnitude of instability X u of symmetric double-peaked solutions of DNLS 
for x < in ID (points) as a function of the strength of nonlinearity for various 
interpeak separations: 5 = 1 (filled circles), S = 2 (filled squares), S = 3 (filled 
diamonds), S = 4 (filled triangles), S = 10 (open circles), S = 15 (open squares), 
S = 20 (open triangles), S 1 = 30 (open diamonds), and S = 50 (crosses). Lines show 
the power-law relation, Eq.(23), derived for relatively large values of |x|. 

This relation is plotted in Fig. 4 for S = 1, 2, 3, 4, 10, 15, 20, and 30 (solid 
lines) along with the corresponding numerical results (points). The agreement 
is very good for \x\ larger than 10 — 20. 

The disappearance of the symmetric branches occurs for non-zero \ (apart 
from the case of S = 1, where the corresponding branch merges with the SP 
branch) when, as in the antisymmetric real eigenvalue splitted from 

the band collides with the pinned eigenvalues at zero. Table 1 shows (fourth 
column) the nonlinearity regime where different branches of symmetric DP 
solutions disappear. It seems that for fixed 5* (S > 1), the antisymmetric 
branches survive until smaller values of \x\- 

As it is shown in the Appendix, close to the anti-continuous limit the stable 
(real) discrete eigenvalue of an antisymmetric DP state has the same depen- 
dence like in Eq. (23) (see Eq. (52)). Numerical simulations confirm that for 
relative large |x| the stable eigenvalues of antisymmetric states have the same 
magnitude with the unstable of the symmetric ones. This result can be used 
for a rough estimate of the nonlinearity value, Xwi where an antisymmetric 
state becomes unstable (through the collision of the discrete eigenvalue with 
the band), by 

M _ 2 = 2 S ' 2 |xM forS>l, or M - 3 = ^2\x~\ for S = 1, (24) 

_ 
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Fig. 5. Frequencies of single- double- and quadruple-peaked stationary solutions of 
DNLS in 2D (points). Dashed lines show analytical expressions obtained for large 
values of |x|. The horizontal line at lo = —4 indicates the lower edge of the band of 
Bloch stationary states, which extends from —4 to 4. The spectrum is antisymmetric 

on x; w(-x) = 

where ^ — 2 — 3) is the lower band edge. The larger the \Xun\ the better 
the estimate, since Eq. (23) is valid for \x\ S> 1. Estimates of Xum resulting 
from the solution of Eq. (24), are shown in the second column of Table 1 inside 
parentheses, next to the numerical results. The relative error is less than 10% 
for S = 1 and S = 2, but it increases for larger S where \Xun\ is getting 
smaller. 



5 Multi-peaked solutions in 2D 

5.1 Frequency spectrum 

Similarly to the ID case, also in 2D there are many families of multi-peaked 
localized stationary states corresponding to discrete levels in the frequency 
spectrum. For large values of \x\ these levels tend to the anti-continuous limit 
spectrum. Fig. 5 shows a part of the frequency spectrum in 2D, which, apart 
from the single-peaked states of lowest frequency, contains many branches 
of double-peaked and quadruple-peaked (QP) stationary solutions at various 
interpeak distances. 

All the calculated DP states (symmetric or antisymmetric, along the lattice 
axes or along the diagonal, and with different interpeak separations) have 
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frequencies around u = |, apart from the symmetric and antisymmetric states 
with their two peaks at neighboring sites (examples are shown at the left 
columns of Fig. 6 and Fig. 7), which have frequencies around uj — | — 1 and 
uj — | + 1, respectively. Therefore, from the branches of the DP solutions of 
Fig. 5, the middle one is highly crowded, tending to the level of Eq.(16) for 
M = 2 at large |x|, while the two external branches are single branches. 

Similar considerations are valid for the branches of quadruplet stationary 
states. What appears as a middle branch of the QP states in Fig. 5 actu- 
ally contains many stationary solutions with frequencies around uj — j . Below 
and above these highly congested branches there are single branches corre- 
sponding to the symmetric (see left column of Fig. 10) and antisymmetric 
(see left column of Fig. 11) solutions, respectively, with their four peaks on 
the corners of the unit cell of the square lattice. 

Multi-peaked stationary states, representative of some of the branches shown 
in Fig. 5, are presented in the following two subsections and their linear stabil- 
ity is discussed. Note that there also exist triple-peaked solutions (with their 
positive or negative peaks along the axes, or along the diagonal, or at random 
sites), which are not discussed here. Only to mention that their correspond- 
ing branches tend to uj = | for large values of \x\, i-e. are in between the 
double-peaked and quadruple-peaked branches shown in Fig. 5. 



5.2 Double- peaked solutions 



DP solutions with the two peaks along a lattice axis (let say the y-axis) are 
presented first. The corresponding branches are calculated starting from the 
initial state 

Vn x ,riy = ~^(^n x ,n 1 ^n y ,n 2 + $n x ,ni$n v ,n2+s) (25) 



for the symmetric and 

) (26) 



for the antisymmetric stationary states, respectively, where S — 1,2, . . . de- 
termines the interpeak separation along the lattice axis. Figs. 6 and 7 show 
examples for S = 1 (left columns), S = 2 (middle columns), and S = 3 (right 
columns). Fig. 5 contains branches of such solutions for S — 1, 2, 3, 4, 5, and 
10. For S = 1 the upper and lower single branches of the DP states of Fig. 5 
are obtained. 
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Fig. 6. 3D plots (first row) and density plots (second row) of double-peaked sym- 
metric solutions (with their peaks along a lattice axis) of DNLS in 2D. Left col- 
umn: interpeak separation S = 1 lattice site, x = —10.5. Middle column: interpeak 
separation S = 2 sites, x = —15. Right column: interpeak separation S = 3 sites, 
X = —14. Cross-sections of the wavefunctions are shown in the third row with points: 
ipn x =n u n y (filled circles), ^ nx=ni+1>riy (filled squares), ip nx = ni +2,n y (filled diamonds), 
V'n ;c =ni+3,n y (filled triangles), where n\ = 20 is the x— coordinate of the two peaks. 
Dashed lines show analytical approximations of the solutions using Eq. (29). 

Branches of DP solutions with their two peaks along the diagonal of the lattice 
axes have been also shown in Fig. 5. These branches are obtained starting from 
the initial state 

^n x ,n y = ~~7^{^n x ,ni^n y ,n2 i $ n x ,^+1$ n y ,n 2 +l) > (27) 



where the plus sign gives the symmetric and the minus the antisymmetric, 
respectively, stationary states. Fig. 5 contains branches of these solutions for 
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Fig. 7. 3D plots (first row) and density plots (second row) of double-peaked an- 
tisymmetric solutions (with their peaks along a lattice axis) of DNLS in 2D. Left 
column: interpeak separation S = 1 lattice site, x = — 14. Middle column: interpeak 
separation S = 2 sites, x = —12. Right column: interpeak separation S = 3 sites, 
X = —14. Cross-sections of the wavefunctions are shown in the third row with points: 
ipn x =n u n y (filled circles), ^ nx=ni+1>riy (filled squares), ipn x = ni +2,n y (filled diamonds), 
V'n ;c =ni+3,n y (filled triangles), where n\ = 20 is the x— coordinate of the two peaks. 
Dashed lines show analytical approximations of the solutions using Eq. (29). 

/ = 1, 2, 3, and 10 and their frequencies are around u = |. Some examples of 
such states are shown in Fig. 8. 

Finally, many other DP solutions exist, with their peaks, of the same or op- 
posite sign, at random lattice sites (not aligned along a lattice axis, or the 
diagonal). These can be obtained starting from the initial state 

^n'n] = "T^OWm'Wna ^ ^n x ,n 1 +sJn v ,n 2 +S y )-, (28) 
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Fig. 8. 3D plots (first row) and density plots (second row) of double-peaked symmet- 
ric and antisymmetric solutions (with their peaks along the diagonal of the lattice 
axes) of DNLS in 2D. Left column: symmetric state with I = 1 and x = —14.1. 
Middle column: antisymmetric state with I = 1 and x = — 11- Right column: sym- 
metric state with I = 2 and x = —13.1. Cross-sections of the wavefunctions are 
shown in the third row with points: tp nx=ni ^ ny (filled circles), ifin x =m+i,n y (filled 
squares), ^=^+2,^ (filled diamonds), ip nx=ni+ 3, ny (filled triangles), ^=711+4,^ 
(open circles, in the second and third column), tp nx=ni+5 ^ ny (open squares, in the 
third column), where n\ = 20 is the x— coordinate of the first peak. Dashed lines 
show analytical approximations of the solutions using Eq. (29). 

with any combination of non-zero integ ers S X ) Sy and S x ^ The cor- 
responding branches are close to the middle branch of the DP solutions of 
Fig. 5. 

As in the ID case, one can obtain approximate analytical expressions for the 
DP solutions of DNLS in 2D, by appropriate superpositions of the single- 
peaked solutions (12). If (rii,n 2 ) is the lattice site of the first peak and {n\ + 
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S x ,ri2 + S y ) of the second one (S x and S y are assumed to be positive), then 
the corresponding approximate solution is 



tb DP 

T Tlx y^y 



1 1 — C 2 

2_ (Anx-ni\+\n y -ri2\ _|_ /_ S g n ^,\S x +S y ^\n x -n 1 -S x \ + \n y -n 2 -S y \ 

^2(1 ±P) 1 + C 2 V 



where P = <,-s 9 n X )^ [£ + ^ = ^ = ™±|^ = (S » = ^ 

1 > 1 6 2 48 

and C = 77T - ; /oX o = o- (31) 

X/2 (x/2) 3 XX 3 



2] /" 'S' :E + 5'-y 



Here also, the individual wavefunctions superimposed in this solution should 
correspond to |, which has been taken into account in the relation (31) pro- 
viding (. The plus and minus signs in Eq. (29) correspond to symmetric and 
antisymmetric states, respectively, apart from the case of positive x an d odd 
S where it is the other way around. The overlap P of the two superimposed 
single-peaked solutions is P = En x ,n y ^}n y lM {\ Wn x }n y l+Sx ' n2+Sy K\ ), where 
ipnxri y ' n2 K f) ^ S ^ ne s °hition in 2D, centered at (ni,ri2) and corresponding 
to |. Cross-sections of the approximate solution (29) are shown with dashed 
lines in the third rows of Figs. 6-8. 

Concerning the stability of the DP stationary states, the picture is similar 
like in ID. For x < 0, or x > and even S, all the symmetric solutions are 
unstable, while the antisymmetric ones are in general (at least for relatively 
large values of |x|) linearly stable. For x > and odd S the reverse is true. 
General arguments are presented in the Appendix showing that close to the 
anti-continuous limit any symmetric (antisymmetric) DP solution should be 
unstable (linearly stable), except when x > and S is odd, where it is linearly 
stable (unstable). Further, this calculation allows to determine the variation 
of the magnitude of instability \ u of the unstable solution with the nonlin- 
earity strength (for relatively large values of |x|), depending on the interpeak 
distance. As it is shown in the Appendix, if the separation of the two peaks 
in the 2D lattice is given by S x and S y , then 

\u = v^T^ 2^)/ 2 Ixl 1 -^, for |X| » I- (32) 



Fig. 9 presents numerical results in the case of x < regarding the magnitude 
of instability of symmetric DP states with different interpeak separations S x , 
S y (points), as well as a comparison with the power-paw (32). 

Besides the pair of stable or unstable discrete eigenvalues, there is also the 
band of eigenvalues extending approximately from ^ — 4 to ^ + 4, except 
for the cases where the two peaks are located in first neighboring sites. Then, 
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Fig. 9. Magnitude of instability X u of symmetric double-peaked solutions of DNLS 
for x < in 2D (points) as a function of the strength of nonlinearity \x\, for various 
interpeak separations. Results for solutions along a lattice axis with S x = 1, S y = 
(filled circles), S x = 2, S y = (filled diamonds), S x = 4, S y = (filled squares), 
S x = 7, S y = (open circles), and along the diagonal with S x = S y = 1 (open 
squares), S x = S y = 2 (filled triangles), and S x = S y = 3 (open triangles), are 
presented. Lines show the power-law relation, Eq.(32), derived for large values of 

Ixl- 

for the symmetric states the band extends from ^ — 3 to -rr + 5, while for 
the antisymmetric states extends from ^ — 5 to ^ + 3. As \x\ decreases, the 
band moves towards zero and the states which are stable at large \x\ become 
unstable when the band collides with the discrete eigenvalues lying on the real 
axis. 

5.3 Quadruple-peaked solutions 

In this subsection QP solutions of high symmetry are presented, where their 
four peaks are on lattice sites forming a square of edge equal to I lattice 
constants. Such solutions may be symmetric, antisymmetric along both lattice 
axes, or symmetric along the one axis and antisymmetric along the other one, 
and at various interpeak distances I. 

The former are calculated from the initial state 

^n x ,n y ~^y^i^n x ,ni^n y ,n2 ~l~ ^n x ,ni^n y ,n 2 +l ~l~ ^n x ,ni+l^n y ,n 2 +l ~l~ ^n x ,ni+l^n y ,n 2 ) • (33) 

Fig. 10 shows such symmetric stationary states for I = 1, 2, and 3. The 
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Fig. 10. 3D plots (first row) and density plots (second row) of quadruple-peaked 
symmetric solutions of DNLS in 2D. Left column: interpeak separation I = 1 lattice 
site, x = —22.5. Middle column: interpeak separation / = 2 sites, x = —32. Right col- 
umn: interpeak separation / = 3 sites, x = —30. Cross-sections of the wavefunction 
are shown in the third row with points: ijj nx=ni ^ ny (filled circles), ij}n x =ni+l,n y (filled 
squares), ^n x = ni +2,n y (filled diamonds), W=ni+3,n y (filled triangles), ipn x = ni +4,n y 
(open circles), i^n x =nx+^,n v (open squares, in the third column), where n\ = 20 is 
the x— coordinate of the first peak. Dashed lines show analytical approximations of 
the solutions using Eq. (36). 

solutions with I = 1 give the single branch with frequencies around uj = j — 2 
in Fig. 5. These solutions are unstable for negative \. Linear stability analysis 
shows three purely imaginary pairs of opposite eigenvalues (two of these pairs 
are degenerate). For the case I — 1 (I > 2), as |%| increases the magnitude of the 
unstable eigenvalues increases (decreases). For I = 2 the unstable eigenvalues 
approach constant values (around ±2.83 the most unstable one and around ±2 
the doubly degenerate eigenvalues) for \x\ 3> 1- However, for relatively large 
X > the symmetric solutions corresponding to odd values of / are linearly 
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Fig. 11. 3D plots (first row) and density plots (second row) of quadruple-peaked an- 
tisymmetric solutions of DNLS in 2D. Left column: interpeak separation I = 1 lattice 
site, x = —38. Middle column: interpeak separation I = 2 sites, \ = —21. Right col- 
umn: interpeak separation I = 3 sites, \ = —26. Cross-sections of the wavefunction 
are shown in the third row with points: ip nx=nijriy (filled circles), ipn x =ni+l,n y (filled 
squares), ip nx =n 1 +2,n y (filled diamonds), rpn x =m+3,n y (filled triangles), ^n x =nx+A,n y 
(open circles), tyn x =n 1 +'o,n y (open squares, in the third column), where n\ = 20 is 
the x— coordinate of the first peak. Dashed lines show analytical approximations of 
the solutions using Eq. (36). 



stable with three real discrete pairs of eigenvalues. The band of eigenvalues 
extends approximately from — 4 to ^ + 4, apart from the case of I = 1, 
where it extends from ^d- _ 2 to ^ + 6. 

The antisymmetric along both axes QP solutions are obtained from alternating 
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signs on neighboring peaks: 



n x ,n y ^f^(dn x ,ni fitly, n 2 ^n x ,ni^n y ,n 2 +l + 0~n x ,ni+l°~n y ,n2+l 0~n x ,ni+l°~n y ,n2) -(34) 



Some examples are shown in Fig. 11. The branch of these solutions with I = 1 
gives the single branch with frequencies around u = j + 2 in Fig. 5. For 
X < these solutions are linearly stable for large values of \x\- In this case 
there exist three discrete real eigenvalues (two of them are degenerate) and 
the band extends from ^ — 4 to ^d. + 4, apart from the case of / = 1, where it 
extends from ^ — 6 to ^ + 2. As long as the discrete eigenvalues are outside 
of the band the solution is linearly stable. When they collide instabilities 
develop. The dependence of the discrete eigenvalues on \x\ is the usual: for 
I — 1 (I > 2) they increase (decrease) with \x\, while for I — 2 they slightly 
vary, approaching constant values at |x| >> 1. Regarding the solutions shown 
in Fig. 11 the last one (for I — 3, in the right column) is linearly stable and the 
other two are unstable. As previously, the opposite sign of x (x > 0) does not 
change anything regarding the stability eigenvalues of symmetric and fully 
antisymmetric QP states, when / is even. On the contrary, for odd I these 
families of solutions interchange stability eigenvalues when x ~^ ~X- 

The last example of QP states, which are symmetric along one lattice axis and 
antisymmetric along the other one, are obtained from the initial state 

^n^n^ ~ ^J-^ i^n x ,ni^n y ,n2 + ^n x ,n\^n y ,n2+l ~ ^n x ,ni+l^n y ,n 2 +l ~~ $n x ,ni+l^n y ,ri2 ) • (35) 

Three cases (for 1 = 1,2, and 3) are presented in Fig. 12. These solutions are 
always unstable since there are two pairs of purely imaginary eigenvalues ±i(3 
and ±ij. There is also a discrete eigenvalue 5 which is real for large values 
of |x|, but when it collides with the band (which extends approximately from 
one more instability is developed. For \x\ 3> 1, ft and 5 tend 
to the same value and they show the typical dependence on \x\' they increase 
(decrease) for I — 1 (I > 2) and they tend to 2 for I — 2. This picture does 
not change for positive x, regardless whether / is even or odd. 

In Fig. 5 branches of solutions of the symmetry obtained from (33) and (34) 
are shown for I — 1, 2, 3, 4, 5, 6, 7, 10, and 11 and from (35) for / = 1, 2, 
3, 4, 5, 10, and 11. Of course there are many more quadruplets having their 
four peaks in a rectangle or in random lattice sites, with any combination of 
signs, which are congested around the middle quasi-degenerate branch of the 
QP solutions of Fig. 5. An example of antisymmetric along both diagonals QP 
solution, forming a square with its edges (of length equal to \^2) not along the 
axes, as the states presented above, but along the diagonals, is reported as a 
quasivortex in Ref. [36] 
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Fig. 12. 3D plots (first row) and density plots (second row) of quadruple-peaked 
symmetric/antisymmetric solutions of DNLS in 2D. Left column: interpeak separa- 
tion I = 1 lattice site, x = —27. Middle column: interpeak separation I = 2 sites, 
X = —27.5. Right column: interpeak separation I = 3 sites, x = —26. Cross-sections 
of the wavefunction are shown in the third row with points: i^n x =n\,n y (filled circles), 
(filled squares), ^ na! =n 1 +2,n s (filled diamonds), ipn x =ni+3,n y (filled trian- 
gles), ipn x=ni +i,ny (open circles), ip nx=ni+5 ^ y (open squares, in the third column), 
where n\ = 20 is the x— coordinate of the first peak. Dashed lines show analytical 
approximations of the solutions using Eq. (37). 



Once more, approximate analytical expressions can be derived for the QP 
solutions in 2D, by superimposing four single-peaked solutions of Eq. (12), 
each one corresponding to f . Regarding the high symmetry solutions presented 
above, if (711,712) is the position of the first peak and I the interpeak distance 
along the lattice axes (/ > 0), then these stationary states can be approximated 
by 
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where P = (-s#nx) ? \ , V 2 (38) 



1 6 4 384 

and C = ^7i-W = ~" (39) 



Plus (minus) signs in Eq. (36) provide the symmetric (fully antisymmet- 
ric along both axes) QP solutions, except when x is positive and I odd, 
where it is the reverse. Eq. (37) gives the solutions which are symmetric 
along one axis and antisymmetric along the other one, like those of Fig. 12. 

P = En x ,n y ^ x }Z' n2 K\)^n x }ny' m+l] (?) is the overla P of two neighboring single- 
peaked states at distance /. The overlap of two single-peaked states across the 
diagonal of the square, En„,n w ^S w 1,n2] (f )^3? +,,n2+,] (! )> ^ equal to P 2 . 
Cross-sections of Eqs. (36) and (37) are shown in the third rows of Figs. 10-12 
with dashed lines. 



6 Conclusions 



Multi-peaked localized excited states are discussed for the one- and two- 
dimensional discrete nonlinear Schrodinger equation. Their numerical calcu- 
lation is achieved by iterations of a simple map, where trivial initial states 
rapidly converge to the desired solution. Examples have been presented of 
symmetric and antisymmetric states with different interpeak separations for 
double-peaked solutions in ID, as well as for double-peaked states along a lat- 
tice axis or along the diagonal and quadruple-peaked states on a square in 2D. 
Analytical approximations and the linear stability of the solutions have been 
discussed. For strong nonlinearities, the symmetric double-peaked states are 
unstable and the antisymmetric linearly stable, except for the case of positive 
nonlinearity and odd interpeak separation, where the situation is reversed. An 
interesting application that such multi-peaked solutions may have, concerns 
their potential use for information encoding and transfer in optical lattices 
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[37]. Multi-peaked solutions have been also discussed in different contexts 
[38,39,40]. 

The classification of these stationary states is based on their origin at the anti- 
continuous limit, even though sometimes their name may be misleading. For 
example, the symmetric double-peaked states in ID and 2D and the symmetric 
quadruplet in 2D with interpeak separations equal to one lattice constant, can 
be viewed as solutions having one peak. However, this classification is very 
useful for organizing the discrete frequency and energy spectrum of DNLS. 
The concept of anti-continuous limit [23] facilitates the interpretation of the 
structure of this spectrum. 

Some of the stationary states presented in this work have been also discussed 
in other studies and different names are attributed to them. The double- 
peaked solution in ID with interpeak separation S = 1 has been named 
Page mode, even mode, or centered-between-sites [41]. Antisymmetric double- 
peaked states in ID with S = 1 or S = 2 are also known as twisted modes 
[42], while their analogues in 2D with their peaks along the diagonal have 
been discussed in [43]. In Ref. [44] the symmetric double-peaked state along 
a lattice axis with S — 1 in 2D and the symmetric quadruple-peaked state 
with interpeak separation / = 1 have been named hybrid mode and Page-like 
mode, respectively. The double-peaked symmetric stationary solution along 
the diagonal with I = 1 in 2D has been also known [45]. 
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Appendix A 

A.l Linear stability of a stationary state 

Considering small deviations 8ip n (t) from a stationary solution of Eq. (5), i.e. 
^n(t) = (jp n + 5ip n (t)) • e~ lU}ot (where ip n is time-independent and real and ujq 
is the corresponding frequency), substituting in DNLS Eq. (1), and linearizing 
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in respect to the complex small perturbations Sip n (t), one obtains 

%d ^T = (_C "° + x ^ S ^ ~ E fyr+s + 2 x ^ 2 n Re(5^ n ), (40) 

where Re(5ip n ) is the real part of Sip n (t). Substituting solutions of the form 
8ipn(t) = a n sin(u;t) + ib n cos(wt) (41) 

in the linearized equation (40), yields the following coupled system 



ua n = (-uo + x4>l)b n ~ b n+s ( 42 ) 
ub n = (-uj + 3x'4''i)a n -J2a n+ s. (43) 

<5 

This implies that the stability eigenvalues u can be obtained as the eigenvalues 
of the 2L x 2L (where L is the total number of lattice sites) matrix M: 




Mi W A 
M 2 \B 




(44) 



where the 2L x 1 column (A, B) T = (a±, . . . , cll, b±, . . . , bi) 7 corresponds to 
the eigenvectors, denotes the L x L zero matrix and Mi, M 2 are the tight- 
binding L x L matrices given through the right-hand-sides of Eqs. (42) and 
(43), respectively. M x and M 2 differ only in their diagonal elements; M lu = 
—u + x'Pf and M 2ii = —ojq + Sxi^f, while the non-diagonal matrix elements 
are zero, except when they correspond to first neighboring sites where they 
are equal to —1. 

From Eq. (41) we see that if there is an eigenvalue of M with non-zero imag- 
inary part, then the stationary solution is unstable. It can be easily verified 
that any stationary solution has a stability eigenvalue uo = 0, with eigenvec- 
tor a n = 0, b n = ip n , since Eq. (43) is trivially satisfied and Eq. (6) provides 
a solution making zero the right- hand- side of (42). The eigenvalue u = is 
doubly degenerate possessing a second eigenvector known as growth mode [8]. 



A. 2 Obtaining the linear stability eigenvalues through a tight-binding Hamil- 
tonian 



Here it is shown that the eigenvalues of the linear stability problem can be 
calculated through the energy eigenvalues and eigenvectors of a tight-binding 
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Hamiltonian 2 with an on-site potential determined by the stationary state 
ip n . In particular, consider the Hamiltonian eigenvalue problem 

HiK = ~ E K+s + U n<t>n = Eu<Pn, with potential U n = x4>l (45) 
s 



where E v and 0^ are the eigenvalues and the corresponding eigenvectors. The 
potential U n has the shape of the stationary state under discussion, multiplied 
by the nonlinearity parameter. Comparing Eqs. (6) and (45) one obtains that 
ip n is an eigenvector 0^ =o of Hi with E u=0 = ujq. 



Expressing the eigenvectors a n and b n of the linear stability system (42), (43) 
in the complete basis <\) v n of the Hamiltonian Hi, i.e. a n = Yliv a v§ v n an d b n = 
J2i/K<t>ni substituting in the system (42), (43), and using Eq. (45) and the 



relation Y.n^n'Pn = <W'> yields that for v ^ is b u = 



E v -Eq 



and 



00 



E,, — Er 



+ 2x 



E(^nC) 2 
n 



h , + 2 X '£a v , E(« 2 « =0.(46) 



These L — l equations (since v ^ 0) provide the 2L — 2 non-zero stability eigen- 
values to = ±\fu~j 2 ~. We see that the eigenvalues of the linear stability analysis 
appear as pairs of opposite sign and the u 2 result from the diagonalization of 
the (L — 1) x (L — 1) matrix 



(E u -E y + 2 X (E u -E )j:(<P n r n ) 2 
+(1 - 5 vy )2 X {E v - Eo)YM 



rCVn i 



S v y (47) 
where v, v' ^ 0. 



K v y is constructed from the eigenspectrum of Hi. Using this matrix, ana- 
lytical results are obtained in the following subsection for the linear stability 
eigenvalues of double-peaked solutions of DNLS for large values of \x\, and 
the Eqs. (23) and (32) are derived. 



A. 3 Application: double-peaked stationary states close to the anti- continuous 
limit 



For a DP stationary solution the potential U n of Eq. (45) is a double well 
for negative x an d a double barrier for positive x- If Ixl ^ 1) then ip n is 

2 a similar procedure has been used in Ref. [16] for the linear stability of a slightly 
more complicated system, viz. the polarons of the adiabatic Holstein model 
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localized at almost two lattice sites, where the two peaks are located, and 
each well (barrier) of U n is very deep (high) and narrow, with a strength 
about |. In this case the Hamiltonian Hi, Eq. (45), is equivalent to a tight- 
binding problem with two equal impurities with large on-site energies |. Then 
the energy spectrum has two discrete eigenvalues (the E Q = u and a second 
one, denoted by E\) and, for an extended system, all the other eigenvalues 
belong to the continuous band from —2c? to 2d [46]. The eigenvectors of the 
continuous spectrum are proportional to ^ and therefore negligibly small at 
any lattice site. Since the sums over the lattice sites n that appear in K v y in 
Eq. (47) contain the tightly localized around two sites wavefunction 0° = ip n , 
they can be neglected when they involve an eigenstate of the continuous. 
The <p v n , <p v n of K v v i do not include the v = and thus there is only one 
localized eigenstate among them; the 0* corresponding to Si. As a result the 
non-diagonal matrix elements of K v y can be neglected, since they contain 
products of two eigenstates 0^0^ and at least one of them it belongs to the 
continuum. The diagonal ones, providing directly the stability eigenvalues uj 2 , 
are 

cj 2 = K UtU = (E v - E ) 2 , for v ^ 1 and (48) 

to 2 = K 1A = (E 1 - E ) 2 + 2 X (E 1 - E ) E(0n0n) 2 , for i/ = 1 (49) 

n 

Eq. (48) gives two bands of real eigenvalues uj = ±\E U — E \, symmetrically po- 
sitioned around zero, and Eq. (49) a discrete pair of eigenvalues uj = ±yK^ 1 . 

The discrete eigenspectrum of Hi is given by 

€ = I = : {g [ :° ] ± (-sgn X ) S 9 [ : o+S] ) ^E ± = 1±£ « 6 ± v T eP, 



where , g^ 0+s ^ represent single-peaked wavefunctions, centered at n and 
n + S, respectively. Here, n, no, and S (S or its components are assumed 
positive) have to be understood as integers in ID and pairs or sums of in- 
tegers in 2D, e.g. n -> (n x ,n y ), n -> (m,n 2 ), S -> (S x ,S y ), {-sgnx) s -> 
(— sgnx) Sx+Sy i etc. The other quantities in the energy eigenvalues E± of Eq. (50) 
are e = (g^\Hi\g^) = {g^^^H^g^ 1 ^^) and the small overlap integrals 

v = (-sgnx) s {gt o] \Hi\gt° +S] ) and p = (-- s 9 n x) S {gt o] \9n 0+S] )- For X nega- 
tive, i.e. attractive potential U n , the upper (lower) signs in Eq. (50) correspond 
to the ground (first excited) state of Hi, while for positive Xi Le - repulsive U n , 
they provide the highest-energy (second highest-energy) state. Note that for 
positive x an d odd S the highest-energy state is the antisymmetric one be- 
cause of the (—sgnx) s term in (50), in contrary to the case of even S. It is 
convenient to distinguish three cases: 
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(i) negative x, 

(ii) positive x an d even S (or even S x + S y ), and 

(iii) positive x an d °dd S (or odd S x + S y ). 

For symmetric DP stationary states ip n , the 0+, -E+ of Eq. (50) correspond to 
(f)° n = *0 n and E'o and the 0~, to 0* and i?i in cases (i) and (ii), while it is 
the other way around in case (iii). For antisymmetric DP states ip n , holds the 
opposite: the 4>~,E_ correspond to 0° and Eq and the (f>+,E + to (j) l n and E\ 
in cases (i) and (ii), and the reverse is valid in case (iii). 

The discrete eigenvalue (49) is dominated by the second term, since E 1 — E 
is a small quantity and x ^> 1. For a symmetric DP solution E\ — E is posi- 
tive in cases (i) and (iii) and negative in case (ii). Therefore, which has 
the same sign as x{E\ — Eq), is negative (positive) in cases (i) and (ii) (in 
case (iii)), meaning that the symmetric DP solution is unstable (linearly sta- 
ble) with a pair of purely imaginary (real) discrete eigenvalues. The situation 
is reverse for an antisymmetric DP stationary state since E\ — Eq has the 
opposite sign now. Thus, when the symmetric DP solution is linearly stable 
(unstable) the antisymmetric one is unstable (linearly stable). These results 
are confirmed by the numerical simulations of sections 4 and 5. The above 
arguments about the linear stability/instability of a symmetric or antisym- 
metric DP solution of DNLS close to the anti-continuous limit, are generally 
applied at any dimension. Therefore, one expects that, for large \x\, also in 
3D an antisymmetric DP stationary state is linearly stable and a symmetric 
unstable, except when x is positive and the interpeak distance S x + S y + S z 
is odd, where they interchange roles regarding their stability. 

The previous discussion can be quantified and using Eqs. (49) and (50) analyt- 
ical expressions are obtained for the discrete pair of eigenvalues in the limit of 
|x| 3> 1. The wavefunctions g n in (50) are given by Eq. (12) for |. Taking into 
account that ( = — ^ one obtains to leading order: P = (—sgnx) s (l + S)( s in 
ID and P = (sgnx) Sx+Sy {l + S X ){1 + S y )( Sx+Sy in 2D [see also Eqs. (19) and 
(30)], the sum J2n(€<Pl) 2 = \, the coupling v = (sgn X ) s (SC 3 ' 1 + x( S ) 
in ID and v = (sgn X ) Sx+Sy [-(Sx + S y + 2S x S y )( s * +s «~ 1 + x( Sx+Sy ] in 
2D, and the on-site energy e = | . Then for symmetric DP states in ID is 

E 1 -E = -2(v - eP) = in case (i), E 1 - E = -2(v - eP) = 

in case (ii), and E\ — Eq = 2{v — eP) = ^1^- in case (iii). The signs are the 

same in 2D, but the magnitudes change to ^ 1+ ^f x y +s y ~i " ■ Therefore Eq. (49) 

yields that in cases (i) and (ii) u 2 = — m ID and u 2 = — ^ 1+ ^$J+$ v -2 ~ 
in 2D, while in case (iii) u 2 is positive with the same magnitude. For the an- 
tisymmetric DP states the E\ — Eq and consequently the u 2 are the same like 
previously, but with opposite sign. As a result the discrete eigenvalues are 

u = ±i2 s / 2 \ X \ l ~ s ' 2 in ID and u = ±ifi + S x S y 2 (^)/2 |x| i-^ in 2D(51) 
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for the (unstable) symmetric DP states in cases (i) and (ii) and the antisym- 
metric in case (iii), while 

uj = ±2 S/2 \ X \ 1 ~ S/2 in ID and uj = ±y/l + S x S y 2( & +^)/ 2 | x | 1 - £ 4 £a i n 2D 

for the (linearly stable) symmetric states in case (iii) and the antisymmetric 
ones in cases (i) and (ii). Eq. (51) provides the instability magnitudes given 
in Eqs. (23) and (32), while Eq. (52) the stable eigenvalue used in Eq. (24). 
By decreasing \x\, as the stationary state ip n and the potential U n of (45) 
become more extended, additional discrete levels of Hi may appear, resulting 
in additional discrete eigenvalues of the linear stability problem. 

For the DP states, E = lu equals, to leading order, to e = |. The next 
corrections are ±(v — eP), which are equal to ± El ~ E " = ± i^fg-i in ID and 
a similar expression in 2D (see above). These corrections are inverse powers 
of Xi apart from the case of S — 1 (or S x + S y = 1). Then, if x < for 
example, the symmetric DP states have E = uo = | + (v — eP) = | — 1 and 
the antisymmetric have E = uoq = | — (v — eP) — | + 1 (see the lower and 
upper branches of DP states in Figs. 1 and 5, while similar is the case for the 
single branches of the QP states presented in 2D). The knowledge of E and 
the continuous band of Hi (from —2d to 2d) provides the bands of stability 
eigenvalues uj = ±\E U — E \, Eq. (48). The latter, for positive uj, extends from 
4r — 2c? to + 2d, apart from the case of S — 1 (or S x + S y — 1), where it 
extends between ^ ± 1 — 2d and ±1 + 2d (the upper signs correspond to 
symmetric DP states in cases i) and ii) and antisymmetric in case iii), while 
the minus signs correspond to the complementary situations). These results 
are in accordance with the numerical observations in sections 4 and 5. 
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